# Name: TransCostRasterDefinition2.py - MCBS 01/29/2015
# Description: 

# import system module
import arcpy
from arcpy import env
from arcpy.sa import *


# Set workspace environment:
# env.workspace = "TransRaster.gdb"
env.workspace = "C:\Users\Marcelo\Documents\Sugarcane\Offline Data\Transportation Raster\TransRaster.gdb"
env.overwriteOutput = True
arcpy.CheckOutExtension("Spatial")

# Precision of raster (in meters):
Precision = 500

# Cost factors for roads:
CRoadDuplicated   = 1
CRoadPavedNDup    = 1.25
CRoadUnpaved      = 1.66
CRoadElse         = 2.5

# Selection criteria for Paved, Unpaved and Waterways
SelectDuplicRoads     = """ TipoPNV_code = 1 """
SelectPavedNDupRoads  = """ TipoPNV_code = 2 OR TipoPNV_code = 7  """
SelectUnpavedRoads    = """ TipoPNV_code = 6 OR TipoPNV_code = 4 """

## Create rasters:
# Duplicated roads
arcpy.MakeFeatureLayer_management("Roads", "DuplicLayer")
arcpy.SelectLayerByAttribute_management("DuplicLayer", "NEW_SELECTION", SelectDuplicRoads)
arcpy.FeatureToRaster_conversion("DuplicLayer", "dummy", "DuplicRaster", Precision)
# Paved not duplicated roads
arcpy.MakeFeatureLayer_management("Roads", "PavedNDupLayer")
arcpy.SelectLayerByAttribute_management("PavedNDupLayer", "NEW_SELECTION", SelectPavedNDupRoads)
arcpy.FeatureToRaster_conversion("PavedNDupLayer", "dummy", "PavedNDupRaster", Precision)
# Unpaved roads
arcpy.MakeFeatureLayer_management("Roads", "UnpavedLayer")
arcpy.SelectLayerByAttribute_management("UnpavedLayer", "NEW_SELECTION", SelectUnpavedRoads)
arcpy.FeatureToRaster_conversion("UnpavedLayer", "dummy", "UnpavedRaster", Precision)

# Raster arithmetic
# Roads
Dup      = Times("DuplicRaster",CRoadDuplicated)
outDup   = Con(IsNull(Dup),CRoadElse, Dup)
Pav      = Times("PavedNDupRaster",CRoadPavedNDup)
outPav   = Con(IsNull(Pav),CRoadElse, Pav)
Un       = Times("UnpavedRaster",CRoadUnpaved)
outUn    = Con(IsNull(Un),CRoadElse, Un)
BestRoad = LowestPosition([outDup, outPav, outUn])
RoadCost = Pick(BestRoad, [outDup, outPav, outUn])

# Save cost rasters:
RoadCost.save("RoadCost1km")

# Destination ports:
#SelectPorts = """ (UF = 'SP' OR UF = 'PR') AND TIPO = 'MARITIMO' """
#SelectPorts = """ (UF = 'SP' OR UF = 'PR' OR UF = 'RJ' OR UF = 'ES') AND TIPO = 'MARITIMO' """
#SelectPorts = """ (UF = 'SP' OR UF = 'PR' OR UF = 'RJ' OR UF = 'ES' OR UF = 'PA' OR UF = 'SC' OR UF = 'RS' ) AND TIPO = 'MARITIMO' """
SelectPorts = """ TIPO = 'MARITIMO' """

arcpy.MakeFeatureLayer_management("PortsPrj", "PortsLayer")
arcpy.SelectLayerByAttribute_management("PortsLayer", "NEW_SELECTION", SelectPorts)
arcpy.FeatureToRaster_conversion("PortsLayer", "UF", "PortsRaster", Precision)

# Execute CostAllocation
# Road:
costAllocOut = CostAllocation("PortsRaster", "RoadCost")
costAllocOut.save("RoadCostAllocation")

# Execute CostDistance
# Road
# Road Cost distance saved as CostDistanceDef3:
outCostDistance = CostDistance("PortsRaster", "RoadCost")
outCostDistance.save("RoadCostDistance")



